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Abstract 

The statistics of low energy states of the 2D Ising spin glass with +1 and -1 bonds are studied 
for L x L square lattices with L < 48, and p = 0.5, where p is the fraction of negative bonds, using 
periodic and/or antiperiodic boundary conditions. The behavior of the density of states near the 
ground state energy is analyzed as a function of L, in order to obtain the low temperature behavior 
of the model. For large finite L there is a range of T in which the heat capacity is proportional 
to T 5 - 33±0 - 12 . The range of T in which this behavior occurs scales slowly to T = as L increases. 
Similar results are found for p = 0.25. Our results indicate that this model probably obeys the 
ordinary hyperscaling relation dv = 2 — a, even though T c = 0. The existence of the subextensive 
behavior is attributed to long-range correlations between zero-energy domain walls, and evidence 
of such correlations is presented. 

PACS numbers: 75.10.Nr, 75.40.Mg, 75.60.Ch, 05.50.+q 



1 



In 1977, Thouless, Anderson and Palmer- (TAP) performed a mean-field theory analysis 
of the ring diagrams which contribute to the free energy of the Ising spin glass.- 1 ^ They 
found that, above the critical temperature T g , the contribution of these ring diagrams was 
subextensive. This means that, while the sum of these diagrams is divergent at T g , their 
contribution at any T > T g can be neglected in the thermodynamic limit.- Therefore, in this 
limit, no signature of the transition is visible in the equilibrium thermodynamic functions 
for T > T g . However, one can still study the critical scaling behavior of finite systems. 

While it is true that hyperscaling is always violated in a mean-field theory, TAP showed 
that a spin glass has severe fluctuations of the order parameter even at the mean- field level. 
Later, it was shown by Sompolinsky and Zippelius^ that the Ising spin glass also violates 
the fluctuation-dissipation theorem. Thus one should not be surprised if it turns out that 
the spin glass does not obey other relations which work for ordinary phase transitions. 

In this work we analyze data obtained from exact calculations of the density of low-energy 
states for finite two-dimensional (2D) lattices. The same data have also been used to study 
the scaling behavior of domain walls for this model.- We will discover that an unusual effect, 
similar to the violation of hyperscaling found in mean-field theory, also occurs in 2D. The 
data were obtained using a slightly modified version of the computer program of Vondrak,- 1 ^ 
which is based on the Pfaffian method. Our data are completely consistent with the data 
of Lukic et al,— ^ which were obtained using the same algorithm. Our analysis of the heat 
capacity is more detailed than theirs, however, and thus we arrive at somewhat different 
conclusions. 

In two dimensions (2D), the spin-glass phase is not stable at finite temperature. Because 
of this, it is necessary to treat cases with continuous distributions of energies (CDE) and 
cases with quantized distributions of energies (QDE) separately.— 1 ^ In this work we will 
study the QDE case. 

The Hamiltonian of the EA model for Ising spins 2 is 

H = -Y,Jij<7i<7j, (1) 

where each spin <7j is a dynamical variable which has two allowed states, +1 and -1. The (ij) 
indicates a sum over nearest neighbors on a simple square lattice of size L x L. We choose 
each bond J^- to be an independent identically distributed quenched random variable, with 
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FIG. 1: (color online) Finite-size scaling of E$/L 2 vs. 1/L 2 . 

the probability distribution 

P{J l3 ) = p5(J t] + 1) + (1 - p)5(J ij - 1) , (2) 

so that we actually set J = 1, as usual. Thus p is the concentration of antiferromagnetic 
bonds, and (1— p) is the concentration of ferromagnetic bonds. Here we will discuss primarily 
the equal mixture case, p = 0.5, but results for p = 0.25 will also be given. 

The ground state (GS) entropy is defined as the natural logarithm of the number of 
ground states. For each sample the GS energy E and GS entropy Sq were calculated for 
the four combinations of periodic and antiperiodic toroidal boundary conditions along each 
of the two axes of the square lattice. When p = 0.5, all four of these types of boundary 
conditions are statistically equivalent. 

Data were obtained for lattices of sizes L = 7, 8, 11, 12, 15, 16, 21, 24, 29, 32, 41 and 48. 
For each L, 500 different random sets of bonds were studied, for each of the four boundary 
conditions. Thus, combining the data for the different boundary conditions, we have 2000 
values of E and So f° r eac h L. 

With the boundary conditions we are using, for which there is no well-defined surface, 
the value of E averaged over samples of the random bonds, is expected to obey 

E /L 2 = e + a e /L 2 (3) 

to lowest order in L. Fig. 1 shows that this works well, and that the value of eo obtained 
from our data is eo = —1.40151 ± 0.00027. In principle, higher order corrections exist,- but 
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FIG. 2: (color online) Finite-size scaling of Sq/L 2 vs. l/L 2 . 

they are not necessary at the level of precision of our data. This agrees with the result found 
by Lukic et alM All statistical error estimates shown in this work represent one standard 
deviation. The best estimate of eo is still the one of Palmer and Adler,— which uses a 
method for which one can go to much larger L, because the entropy is not calculated. 

The finite-size scaling behavior of So is slightly more complex. Lukic et alr& used a single 
correction-to-scaling term, with an exponent — (2 + 0s). From a fundamental viewpoint,— 
however, when 05 is positive the natural form to use when adding another fitting parameter 
is 

S /L 2 = s + a s /L 2 + b s /L 4 . (4) 

In Fig. 2 we see that this form works well, and gives a value of Sq = 0.07211 ± 0.00015. 
This value is slightly higher than the one quoted by Lukic et al, but the difference comes 
primarily from the different form of the fitting function rather than from differences in the 
data. By comparing with the work of Bouchaud, Krzakala and Martin,— one sees that Lukic 
et al. have made a sign error, and that their fit actually uses a negative value for 0^, which 
is incorrect.— 1 ^ 

While our values of the energy and entropy of the GS of finite LxL lattices for p = 0.5 are 
generally consistent with those of other workers, our results for L = 32 differ substantially 
with those reported by Blackman and Poulter.— (See Figs. 7 and 8 of their paper.) The 
origin of this discrepancy is unclear, but it appears to be too large to be explained by the 
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FIG. 3: (color online) Scaling of 5*1 — So with L. 

different boundary conditions used by them. Their numbers of samples computed are rather 
small, and it may be that they have simply underestimated their statistical errors. However, 
their algorithm, unlike the one used here, does not use exact integer arithmetic to calculate 
the partition function. Therefore, it is likely that they have a problem with roundoff errors. 
In a strongly correlated system such as the one we are studying, substantial roundoff errors 
can result in distributions which are too narrow. 

In order to obtain information about the low temperature behavior, it is useful to study 
the scaling with L of 5i — So, which is the logarithm of the ratio of the degeneracies of the 
lowest excited state and the GS.— We found that 

av(ln(5i - 5b)) = pln(ln((L 2 )!)) + 0.528 ±0.011, (5) 

with 

p = 0.1921 ±0.0015 (6) 

gives an excellent fit for L > 10, as shown in Fig. 3. av() is a configuration average over 
random samples. The points for L = 7 and 8 (not shown in the figure) are below the fitted 
line, due to corrections to scaling at small L. 

The choice ln(ln((L 2 )!)) may appear arbitrary to the reader, but it was suggested by 
the behavior of the fully frustrated 2D Ising model.— In principle, if one could go to 
very large values of L, one could obtain p by plotting the data against 21n(L). From 
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FIG. 4: Histogram of the distribution of Si — Sq for L = 48. 



Stirling's approximation one sees immediately that the difference between using 2 ln(L) and 
ln(ln((L 2 )!)) is a logarithmic correction to scaling. This logarithmic correction appears to be 
present in the data, however, and a much better fit is obtained if one does things as shown 
here. 

If one uses 2 ln(L — 3) an excellent fit over the range of the data is obtained. However, 
this seems completely artificial to the author. In any case the value of p = 0.1948 ± 0.0008 
which one finds from this form is close to the one shown in Fig. 3. (The reason why the 
statistical error in this number is so small is that no contribution from the uncertainty in 
the fitting parameter "3" is included.) 

The reason for taking the configuration average of In (Si — So) rather than taking the 
logarithm of av(Si — So) is that in this way we find the most probable value.— £L The 
probability distributions for Si — S are highly skewed, and the most probable value scales 
differently with L than av(Si — So) does. To illustrate this point, in Fig. 4 we show a 
histogram of the distribution of Si — So for the L = 48 lattices. If one plots the data using 
av(Si — So), one finds an apparent value for p of 0.233(3). Using the median value gives 
0.222(3). It is the typical or most probable value which is the experimentally observable 
quantity, as established by Edwards and Anderson^ for the spin glass. 

From this analysis, we obtain the typical value of Si — S to be 



Si -S = f(L) 



A[\n((L 




(7) 
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with A = 1.696 ± 0.019, or, using Stirling's approximation, 

f(L) « A[L\2\n{L) - l)r . (8) 

It follows immediately that the scaling of Si — So with L is approximately a power law, 
with an exponent close to 0.4, times ln(L). This variation with L is much more rapid than 
the hypothesis of Wang and Swendsen,— who argued for a dependence like 41n(L). To this 
extent, it agrees with the claims of Jorg et al.— 

To obtain the actual behavior of the low temperature specific heat, we must carry the 
analysis further. The heat capacity of a sample of size L x L at temperature T is given by 

C(L,T) = ((E(L)-(E(L))) 2 )/T 2 , (9) 

where the angle brackets indicate a thermal average, and we are using units in which Boltz- 
mann's constant is 1. 

Writing the partition function of a finite sample with periodic boundary conditions ex- 
plicitly gives 

l-Eol/2 

Z(T) = e MSn -So - An/T) . (10) 

ra=0 

The heat capacity is then 

\E \/2 

C(L,T) = (T 2 Zy x J2 16(n - n,) 2 exp(S n - S - An/T), (11) 

n=0 

where n* is the value of n for which the argument of the exponential has its maximum for 
a given sample at temperature T. 

The average values of S n — So for small values of n are shown in Fig. 5, over our full range 
of L. The slope defined by these points, omitting the n = 1 points, is plotted versus 1/L in 
Fig. 6. The limiting value of this slope for large L obtained from this plot is found to be 

ip = 0.842 ±0.003. (12) 

This means that for n <^ L 2 

S n -S w f{L)n*, (13) 

and implies 

ML.T)) 1 -* » i>f(L)T/ 4 . (14) 
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FIG. 5: (color online) Scaling of av(S n — So) with L, for small values of n. The axes are scaled 
logarithmically. 
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FIG. 6: (color online) Slope of av(S n — So) vs. 1/L, for n = 2 to 8. 

This can only be valid, however, if < <C L 2 . If we take the limit T — > 0, holding L 
fixed, then — > 0. Thus the limiting low temperature behavior of C(L,T), for any fixed L 
is proportional to exp(— 4/T), as it must be. We expect to see this behavior when T < T\, 
where 

T X {L) = 4[L 2 (21n(L)-l)]-7(iM) (15) 
is the temperature where n* = 1. We have found a positive value for p, which means that 
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T\ — > as L — > oo. 

The reason for omitting the n — 1 points shown in Fig. 5 from the fits is that they all lie 
well below the straight lines. The quantity Si — Sq does not behave in the same way that 
the other S n+ i — S n do. The author understands this effect by analogy with the well-known 
behavior of random matrices. The gap at the band edge is special, because it only feels level 
repulsion from one side. 

Substituting our expressions for n* and S n — So into Eqn. (10) gives 

C(L,T) = 16T , (16) 

where 

g(n,L) = n i '[A{L 2 {2ln{L) -l)) p - 4n 1_ ^/T] . (17) 

When we try to take the limit L —>■ oo holding T fixed, we get a surprise. The exponent 
p/(l — iji) is 1.216 ± 0.033. Because this exponent is greater than than 1, the power-law 
behavior described by the exponent of Eqn. (12) is only valid for T < T x , where T x must 
go to zero as L increases, n* cannot become larger than L 2 \ This condition requires that, 
when L — > oo, T x must also go to zero at least as fast as 

T X (L) ~ 4L" 2 ^- 1 )(21n(L) - l)-"/{^A). (18) 

Although we do not have data to show that T x actually behaves precisely in this way, it is 
at least plausible that T x goes to zero more slowly than T\ as L increases, since ip < 1. 

What this means is that the singularity we are studying is subextensive, just as the 
thermal singularity above T g is in the TAP mean-field theory.-^ It also means that for L 
large, but finite, we expect there exists a temperature regime T\ T T x in which the 
scaling behavior controlled by this singularity is observable. 

p controls the thermal behavior in the temperature range < T < T 1; and ip controls the 
behavior in the range T\ < T < T x . Therefore, these exponents are independent. A simple 
scaling relation between exponents defined in different ranges of T which have independent 
behaviors is impossible. This statement is not in contradiction with the fact that the value 
of T x clearly depends on both p and ip. The entire procedure used here is quite similar to 
the theory of nested boundary layers. 22 

Since (E(L)) is essentially 4n*(L,T), the heat capacity for T\ < T < T x is easily seen to 
be proportional to T^ 1 "^, which is T 5 - 33±a12 . Because T x as L — > oo, this behavior 



disappears in the thermodynamic limit. The exponent 2(p + if> — 1) is 0.068 ± 0.009. This 
is small, so T x is going to zero quite slowly. Thus the power-law behavior of C(L, T) should 
be visible for macroscopic values of L. Note that this effect is not caused by our choice of 
logarithmic averaging of Si — So, since the use of simple averaging would give a larger value 
for p. 

Although our statistical errors are small, the estimate of p depends on our choice of the 
finite-size scaling fitting function. Notice that the estimate of the scaling exponent a for 
the T dependence of C(L,T) depends only on if), and is independent of p. Therefore, our 
estimate 

a = -5.33 ±0.12 (19) 

is independent of whether T x — > as L — > oo. 

All of the calculations for p = 0.5 described above were repeated for p = 0.25. Using the 
same procedures as discussed above, we find for p = 0.25 the exponents p = 0.1874 ± 0.0019 
and i/j = 0.8527 ± 0.0017. Therefore we obtain 2(p + ip - 1) = 0.080 ± 0.007 and a = 
—5.79 ± 0.08. These results are quite consistent with universality of the critical exponents, 
since the quoted statistical errors do not include any allowance for errors in the assumed 
scaling forms. 

Recently, Jorg et ai— have claimed that a power-law behavior of C(L, T) is evidence that 
the QDE is in the same universality class as the CDE. However, they have not calculated a 
directly. They have calculated the correlation length exponent v « 3.5, and assumed that 
a could be obtained via the modified hyperscaling relation of Baker and Bonner M The fact 
that our value of a is not close to —7 shows that this relation is not obeyed. Our value seems 
to indicate that the ordinary hyperscaling relation, dv = 2 — a, is obeyed, a has never been 
calculated directly for the CDE, so we cannot say whether the values of a are the same for 
the QDE and the CDE. 

Finally, we discuss the origin of the subextensive singularity. Such behavior in a 2D model 
probably requires the existence of some kind of long-range interactions. Such interactions 
are not present explicitly in our Hamiltonian, Eqn. (1), but they may arise spontaneously. 
Since domain walls are extended objects, it would not be very surprising for interactions 
between domain walls to have long range, especially at T = 0. 

Using the same computer program which was used here to obtain the heat capacity and 
additional procedures described in a recent publication,- we have calculated the average 
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FIG. 7: (color online) Scaling of average entropy of zero-energy domain walls for lattices of size 
L x M vs. M/L 125 . These domain walls run across the lattice in the short direction, which has 
length L. 

domain-wall entropy for zero-energy domain walls on lattices of size L x M, where L < M. 
Remarkably, the average domain-wall entropy for the zero-energy domain walls which run 
across the lattice in the short (L) direction scales to zero exponentially in the variable 
M/L 1 ' 25 . This is shown in Fig. 7. The exponent 1.25 is suggestive of the relation recently 
proposed by Amoruso, Hartmann, Hastings and Moore,— which gives a value of 1.25 for 
the fractal dimension of domain walls for this model. From the data displayed here we can 
say that this exponent must be 1.25 ± 0.05. Because their entropy scales to zero so rapidly, 
these zero-energy domain walls must be highly correlated. 

This effect is strong evidence for long-range interactions between the zero-energy domain 
walls. It does not occur for domain walls of other energies. Amoruso et al. do not explicitly 
specify that the behavior of the zero-energy domain walls should be special. However, this 
was suggested by the work of Wang, Harrington and Preskill.— This domain-wall entropy 
calculation will be described more fully in a subsequent publication.— 

In this work we have calculated in detail the low temperature thermal behavior of the 
2D Ising spin glass with an equal mixture of +1 and —1 bonds. We have found that this 
behavior bears a strong qualitative resemblance to the behavior found in the TAP mean- 
field-theory analysis. For finite L there is a range of T for which C(L, T) is proportional to 
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T 5 33 . However, this behavior disappears slowly as L — > oo. This subextensive behavior is 
attributed to correlations between zero-energy domain walls. 
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